clear; clc; close all;

% 参数范围和步长
P_vals = linspace(-5, 5, 100);
Q_vals = linspace(-5, 5, 100);

% 固定参数
q = 1; % 整数阶
a = 1;
b = 1;
x0 = 2; y0 = -2;
N_total = 300; % 轨迹长度
N_transient = 10  0; % 丢弃暂态

% 近似熵参数
m = 2;
r_factor = 0.2;

FraAppEn_map = zeros(length(Q_vals), length(P_vals));

parfor iP = 1:length(P_vals)
    P_val = P_vals(iP);
    temp_row = zeros(1, length(Q_vals));
    for iQ = 1:length(Q_vals)
        Q_val = Q_vals(iQ);

        % 计算轨迹
        [~, y] = SCLMM(q, a, b, P_val, Q_val, x0, y0, N_total);

        % 丢弃暂态，取后半段 y 序列
        y_data = y(N_transient+1:end);

        % 计算容忍度 r
        r = r_factor * std(y_data);

        % 计算分数阶近似熵
        temp_row(iQ) = fractionalApproxEntropy(y_data, m, r);
    end
    FraAppEn_map(:, iP) = temp_row';
end

% 绘制彩色图
figure;
imagesc(P_vals, Q_vals, FraAppEn_map);
axis xy;
colormap(jet);
colorbar;
xlabel('P');
ylabel('Q');
title(sprintf('Fig.14(a) 2D 彩色 FraAppEn 图 (q=%.2f, a=b=1)', q));
set(gca, 'FontSize', 12);

% ------- SCLMM函数 -------
function [x, y] = SCLMM(q, a, b, P, Q, x0, y0, N)
    x = zeros(1, N);
    y = zeros(1, N);
    x(1) = x0; y(1) = y0;

    w = zeros(1, N);
    for k = 1:N
        w(k) = exp(gammaln(k + q - 1) - gammaln(k));
    end

    delta_x = zeros(1, N);
    delta_y = zeros(1, N);

    for n = 2:N
        delta_x(n-1) = a * sin(P * y(n-1)) * sin(P * x(n-1)) - x(n-1);
        delta_y(n-1) = b * cos(Q * x(n-1)) * cos(Q * y(n-1)) - y(n-1);

        sum_x = 0; sum_y = 0;
        for i = 1:n-1
            sum_x = sum_x + w(n - i) * delta_x(i);
            sum_y = sum_y + w(n - i) * delta_y(i);
        end

        x(n) = x0 + sum_x / gamma(q);
        y(n) = y0 + sum_y / gamma(q);
    end
end

% ------- fractionalApproxEntropy函数 -------
function apEn = fractionalApproxEntropy(data, m, r)
    N = length(data);
    phi = zeros(1, 2);

    for k = m:m+1
        C = zeros(1, N-k+1);
        for i = 1:N-k+1
            temp = data(i:i+k-1);
            count = 0;
            for j = 1:N-k+1
                if max(abs(temp - data(j:j+k-1))) <= r
                    count = count + 1;
                end
            end
            C(i) = count / (N-k+1);
        end
        phi(k - m + 1) = (1/(N-k+1)) * sum(log(C));
    end

    apEn = abs(phi(1) - phi(2));
end

